restore,'grid.sav'

default, ivtk, 0
default, ivdf, 0

nr=(size(r))[1]
nth=(size(th))[1]
nphi=(size(ph))[1]
file='fort.11'
n=220D
uu=fltarr(nphi,nth,nr)
vv=fltarr(nphi,nth,nr)
ww=fltarr(nphi,nth,nr)
tt=fltarr(nphi,nth,nr)
pp=fltarr(nphi,nth,nr)
bxx=fltarr(nphi,nth,nr)
byy=fltarr(nphi,nth,nr)
bzz=fltarr(nphi,nth,nr)

w=fltarr(nphi,nth,n)
bx=fltarr(nphi,nth,n)
ns=0L

openr,lun,file,/get_lun,/f77_unformatted
for ns=0,n-1 do begin 
  ; while (not eof(lun)) do begin
  readu,lun,uu
  readu,lun,vv
  readu,lun,ww
  readu,lun,tt
  readu,lun,pp
  readu,lun,bxx
  readu,lun,byy
  readu,lun,bzz

  w[*,*,ns]=reform(ww[*,*,44])
  bx[*,*,ns]=reform(bxx[*,*,43])
  print,'reading data from file: ',file, ns
  ;   ns++
  ; endwhile
endfor
free_lun,lun

if (ivtk ne 0) then begin
 rectilinear_grid=1
; Open the VTK file for write.
 file = 'work.vtk'
; Find the dimensions of the data.
 x=ph
 y=th
 z=r*1000.
 dx=x[1]-x[0]
 dy=y[1]-y[0]
 dz=z[1]-z[0]
 data_type = strlowcase(size(x, /tname)) 
 print,data_type
 nx=n_elements(x)
 ny=n_elements(y)
 nz=n_elements(z)
 ntot=nx * ny * nz
 dimensions= [nx, ny, nz]
 origin=[x[0], y[0], z[0]]
 spacing=[dx,dy,dz]
 ndim=n_elements(where(dimensions gt 1))
	       
; Open the VTK file for write.
 print, 'Writing ', strtrim(file), '...'
; openw, lun, file, /get_lun, /swap_if_little_endian
 openw, lun, file, /get_lun
		   
; Write the header information.
 printf, lun, '# vtk DataFile Version 2.0'
 printf, lun, 'EULAG Data'
 printf, lun, 'ASCII'
 if (rectilinear_grid eq 1) then begin
  printf, lun, 'DATASET RECTILINEAR_GRID'
  printf, lun, 'DIMENSIONS ', dimensions
  printf, lun, 'X_COORDINATES ', nx, ' ', data_type
  printf, lun, x
  printf, lun, 'Y_COORDINATES ', ny, ' ', data_type
  printf, lun, y
  printf, lun, 'Z_COORDINATES ', nz, ' ', data_type
  printf, lun, z
 endif else begin
  print, 'writing header'
  printf, lun, 'DATASET STRUCTURED_GRID'
  printf, lun, 'DIMENSIONS ', dimensions
  printf, lun, 'POINTS ', ntot
  printf, lun, x
  printf, lun, y
  printf, lun, z 
 endelse

 printf, lun, 'POINT_DATA ', ntot
; Write out each data field.
; Write scalar field.
 print, 'SCALARS ', strlowcase('bx') , ' ', data_type, '...'
 printf, lun, 'SCALARS ', strlowcase('bx') , ' ', data_type
 printf, lun, 'LOOKUP_TABLE default'
 printf, lun, bxx

 print, 'SCALARS ', strlowcase('by') , ' ', data_type, '...'
 printf, lun, 'SCALARS ', strlowcase('by') , ' ', data_type
 printf, lun, 'LOOKUP_TABLE default'
 printf, lun, byy
 
 print, 'SCALARS ', strlowcase('bz') , ' ', data_type, '...'
 printf, lun, 'SCALARS ', strlowcase('bz') , ' ', data_type
 printf, lun, 'LOOKUP_TABLE default'
 printf, lun, bzz
 
 ; print, 'VECTORS ', strlowcase('bfield'), ' ', data_type
; printf, lun, 'VECTORS ', strlowcase('bfield'), ' ', data_type
;; printf, lun, swap_endian(bxx,byy,bzz, /swap_if_big_endian)                 
; printf, lun, bxx,byy,bzz                 


 close, lun
 free_lun, lun                     
endif

; VAPOR DATA COLLECTION
;
if (ivdf eq 1) then begin
 isph=1
 default, tinit, '0'
 default, tend, '1'
 default,variables,['u','v','w','bx','by','by']
 if (isph eq 1) then begin 
  varnames=['vphi','vtheta','vr','bphi','btheta','br']
 endif else begin
  varnames=['vx','vy','vz','bx','by','bz']
 endelse

 cd, current=c
 readpath=c
 spawn,'mkdir VAPOR'
 writepath=c+'/VAPOR'

; Find the dimensions of the data.
 x=ph
 y=th
 z=r*1000.
 dx=x[1]-x[0]
 dy=y[1]-y[0]
 dz=z[1]-z[0]
 data_type = strlowcase(size(x, /tname)) 
 nx=n_elements(x)
 ny=n_elements(y)
 nz=n_elements(z)
 ntot=nx * ny * nz
 origin=[x[0], y[0], z[0]]
 spacing=[dx,dy,dz]

 if (isph eq 1) then dim = [nz,ny,nx] else dim = [nx,ny,nz] 

 ; coarsened approximations to create
 num_levels = 1   
 ; create a new VDF metadata object
 mfd = vdf_create(dim,num_levels)
 ; number of timesteps in the dataset
 timesteps = tend-tinit+1                                                                                     
 vdf_setnumtimesteps, mfd,timesteps
 ; set the names of the variables the data set 
 numvar = n_elements(varnames)
 print,'number of variables to write =', numvar
;
 vdf_setvarnames, mfd, varnames
 ; set extents of the volume
 if (isph eq 1) then begin 
   vdf_setcoordtype,mfd,'spherical'
   ; order for spherical coordinates phi(z),theta(y),r(x) 
   extents = [-180+(180.*x[0]/!pi),   -90+(180.*y[0]/!pi),    z[0], $ 
              -180+(180.*x[nx-1]/!pi),-90+(180.*y[ny-1]/!pi), z[nz-1]]
   print,extents
   ;new variables
   vphi=fltarr(nx,ny,nz)
   vtheta=fltarr(nx,ny,nz)
   vr=fltarr(nx,ny,nz)
   br = fltarr(nx,ny,nz)
   bt = fltarr(nx,ny,nz)
   bp = fltarr(nx,ny,nz)
   last_dim = nz-1
  endif else begin
   extents = [x[0], y[0], z[0], x[nx-1], y[ny-1], z[nz-1]]
   last_dim = nz-1
  endelse
 ;
 vdf_setextents, mfd, extents
 ; Store the metadata object in a file for subsequent use
 vdffile = 'work.vdf'
 vdf_write, mfd, vdffile
 ; Destroy the metadata object. We're done with it.
 vdf_destroy, mfd

 ; Writing data
 ;
 itt = 0                           
 for it = tinit,tend do begin ; time loop
   for nvar = 0, numvar-1 do begin ; variables loop
    dfd = vdc_bufwritecreate(vdffile)             
    vdc_openvarwrite, dfd, itt, varnames(nvar), -1
    ; Write (transform) the volume to the data set one slice at a time
    if (keyword_set(sph)) then begin
     case nvar of      
      1: tmp_var = vphi       
      2: tmp_var = vtheta
      3: tmp_var = vr    
      4: tmp_var = bphi       
      5: tmp_var = btheta       
      6: tmp_var = br   
     endcase           
    endif else begin     
     case nvar of     
      1: tmp_var = vx       
      2: tmp_var = vy
      3: tmp_var = vz    
      4: tmp_var = bx       
      5: tmp_var = by       
      6: tmp_var = bz   
     endcase           
    endelse              
    for z = 0, last_dim do begin                                                                             
     vdc_bufwriteslice, dfd, float(tmp_var[*,*,z])
    endfor
    ; close variable / timestep
    vdc_closevar, dfd
   endfor
   itt++
 endfor
 ;       Destroy the "buffered write" data transformation object. 
 ;       We're done with it.
 cd, readpath
 vdc_bufwritedestroy,dfd
endif




end
